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Abstract 

We adapt the VMR (valleys-to-mountains reflections) algorithm, orig- 
inally devised by us for simulations of SOS models, to the BCSOS model. 
It is the first time that a cluster algorithm is used for a model with con- 
straints. The performance of this new algorithm is studied in detail in 
both phases of the model, including a finite size scaling analysis of the 
autocorrelations. 



PACS numbers: 68.35.-p, 68.35.Rh, 64.60.Fr, 64.60.Ht, 05.50.+q 



Introduction 

Cluster algorithms are one of the most promising new ideas to overcome the problem 
of critical slowing down (CSD) in computer simulations of statistical physics models 
[1], Qj. Each time they were sucessful, the dramatic improvement in performance 
came by integrating into the algorithm the relevant physical properties of the model 
under study In the case of solid-on-solid (SOS) models J3|, we devised a new cluster 
algorithm Q that is based on reflecting "valleys" and "mountains" of the fluctuating 
SOS surface. We shall call this algorithm the VMR algorithm (valleys-to-mountains- 
reflections) . 

Up to now, no cluster algorithms were devised for models with constraints on the 
configurations. In this letter we present a first example where this is achieved: the 
adaptation of the VMR algorithm to the BCSOS model. 

The BCSOS model is defined as follows [|[ |J. On each site i of a square lattice 
there is an integer- valued variable hi, which can be thought of as the height of a sur- 
face above the two-dimensional lattice. Nearest neighbor variables are constrained 
to a height difference of ±1. Next to nearest neighbors (i.e. diagonal neighbors) 
can thus have a height difference of ±2 or 0, and the ±2 is assigned the smaller 
Boltzmann weight. The partition function is then 



where are pairs of diagonal neighbors (J2' is the constrained sum). The model 
undergoes a roughening transition at K r = In 2. At large K the SOS surface is 
smooth, at small K it is rough. The BCSOS model is equivalent to the F-model 
@j Hi 01 j which is a particular case of the 6- vertex model (the coupling K in eq. flU) 
is chosen such as to conform with standard 6- vertex notations). 

The idea for the VMR algorithm came from the picture of the SOS configurations 
as landscapes with mountains and valleys. We first choose a horizontal reflection 
plane. The clusters are, roughly, the connected regions above the plane (mountains) 
and below it (valleys). Large scale changes of a configuration are performed by "flip- 
ping" valleys or mountains, i.e. by reflecting them through the plane independently 
and with an appropriate probability. Then a new reflection plane is chosen, the 
reflection procedure repeated, and so on. 

In one way or another the reflection idea had turned up in other cluster algo- 
rithms H |9|, |2|. In the case of the VMR algorithm for SOS models however, the 
crucial ingredient in eliminating CSD altogether turned out to be not only the re- 
flection, but the precise procedure for choosing the reflection plane ||. It is at this 
stage that the particular features of the SOS models were incorporated into the 
algorithm. 

We shall see that the VMR algorithm for the BCSOS model eliminates CSD only 
half-way. It is nevertheless of considerable algorithmic significance to achieve even 
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that much for a model with constraints. Besides, the improvement in performance 
was enough for the purpose of a high precision determination of the roughening 
transition in several SOS models |TIJ. In that study the renormalization group flows 



of the SOS models were matched to that of the critical trajectory of the BCSOS 
model. 



The VMR algorithm for the BCSOS model 

We shall describe a "single cluster flip" algorithm || |2j. A step consists in growing 
a cluster around a seed, then flipping it. The seed is a randomly chosen site Iq. 

For a given seed io, the first thing is to choose an integer reflection plane M. 
We tried two choices of M. One is to set M = h io ± 1, the +1 or —1 being 
chosen randomly with probability one half. The other is to choose another random 
site jo, and set M = hj . Both choices gave the same dynamical exponent z (see 
definition below). The smallest autocorrelation times were obtained with the second 
procedure, in the case that jo is taken on the sublattice that does not contain i (we 
have an even and an odd sublattice; if on the even sublattice all heights are even 
integers, then on the odd sublattice they are odd integers, and vice versa). 

Starting from the seed, which is the first site in the cluster, we freeze and delete 
links with probabilities to be given below. We have to consider both nn links (nearest 
neighbor) and nnn links (next to nearest neighbor). A new site is taken into the 
cluster if it is connected by a frozen link to a site already in the cluster. We repeat 
the freeze/delete procedure until all links pointing from sites in the cluster either 
end in the cluster or are deleted. After the cluster is finished, we reflect it (flip it) 
through the reflection plane M: 

hi^2M- hi (2) 

for all sites % in the cluster. 

The freeze/delete probabilities are defined as follows (for each link the sum of 
the freeze and delete probability is one). For the nn link < i,j > we have two 
situations. If either hi = M or hj = M, the link is always deleted. Otherwise (i.e. 
if both sites are on the same side of the reflection plane) the link is always frozen. 
Thus it is ensured that the constraint on the nearest neighbors is not violated by 
the reflection. 

As a consequence of the freeze/delete procedure for nn links, the nnn link [i,j] 
is automatically frozen in all but three situations (remember that \h{ — hj\ is zero 
or two). First, if hi and hj are on different sides of the reflection plane, the nnn link 
[i,j] is always deleted. Second, if either hi = M or hj = M, [i,j] is also deleted 
with probability one. Third, if hi = hj = M ± 1, the nnn link [i,j] is deleted with 
probability exp(— K). In eq. ([[]) we wrote the nnn interaction in the form of the 
discrete Gaussian model. It can be easily seen that the delete probabilities used for 
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the nnn links are indeed the same as those defined in [|J for the discrete Gaussian 
model without constraints. 

Notice that for the case of the BCSOS model, the clusters are precisely valleys 
or mountains cut by the reflection plane from the SOS surface. In the case of the 
discrete Gaussian model [f§], this was only approximately so. 



Autocorrelations, sublattice effects 

We shall call "autocorrelation time" r the quantity often referred to in the litera- 
ture || as the "exponential autocorrelation time". At large "Monte Carlo times" 
t the autocorrelation function of some physical quantity decays exponentially as 
exp(— t/r), which is related to the slowest mode in the Markov process used for the 
simulation. 

In practice, for some physical quantities the autocorrelations immediately en- 
ter the exp(— t/r) regime, while for other quantities this only happens after the 
autocorrelations have decayed much faster for a while. Of course, the theoretical 
expectation is that at very large t the autocorrelations decay with the same r for 
all quantities (except if prevented by some conservation law). 

For the BCSOS model we found that quantities defined on one sublattice tended 
to behave as exp(— t/r) already for t of the order of a few to a few tens of cluster 
flips. On the other hand, averages over the two sublattices tended to enter this 
"clean exponential" regime only after the autocorrelation functions had decayed 
over several orders of magnitude. 

As an example we discuss here the energy density, which is, up to a constant 
factor, < V" 1 J2[i,j](hi — hj) 2 >, with V the number of points on one sublattice, and 
the sum over nnn pairs in that sublattice. The values of the energy density 
on the two sublattices are so strongly anticorrelated that the statistical error of the 
average over the two sublattices is a whole order of magnitude less than the error of 
the one-sublattice quantity. 

These strong anticorrelations can be understood by regarding the four points at 
the corners of an elementary plaquette of the lattice. In the case that the height 
difference of the nnn variables at one pair of diagonally opposite sites is two, the 
constraint on the nn pairs implies that for the other pair of diagonally opposite sites 
the height difference is necessarily zero. The anticorrelations stem from the fact 
that the two nnn pairs discussed here lie on different sublattices. 

In what follows, we shall present the results for the autocorrelation time of 
the one-sublattice energy. At each value of K and for each lattice size quoted, 
the statistics was between 250000 and 800000 clusters (vectorization of the cluster 
algorithm fill] was °f g rea t help). We used periodic boundary conditions in the 



directions of the coordinate axes of the sublattices (these axes are rotated by 45° 
with respect to the coordinate axes of the original lattice). Each sublattice has a 
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size of L 2 , so our lattice has volume 2L 2 . For our study we took L to be a power of 
2, and 8 < L < 256 (in the rough phase the largest L was 128). 

As opposed to the one-sublattice energy, for the energy averaged over both sub- 
lattices the exp(— t/r) regime is not reached before the autocorrelations become zero 
within errorbars, if the statistics is a few hundred thousand of clusters. In some runs 
with ten times as many clusters we checked that the energy averaged over the two 
sublattices does indeed have the same r. 

The values for r we shall quote are in units of sweeps, i.e. we compute r in units 
of clusters first, then multiply it by the average cluster size and divide it by 2L 2 . 
One unit then amounts to a work proportional to the volume. As opposed to other 
cluster algorithms as e.g. the single cluster algorithm for the Ising model || [| the 
average cluster size is quite stable. If we choose the reflection plane to be the height 
at a random site on the sublattice not containing the seed, the average cluster size 
is around 0.35 of the volume. 

Autocorrelation study in the rough phase and at K r 

Since the SOS surface thickness increases only logarithmically with the lattice size, 
the SOS surface has in reality shallow valleys and small mountains. As the temper- 
ature increases (K decreases) the valleys become deeper and the mountains higher. 
This leads us to the expectation that the higher the temperature, the better our 
algorithm should perform. 

We studied the autocorrelations at the roughening transition K r = In 2 and at 
the point K = K r /2 which is deep in the rough phase. Since in the rough phase the 
correlation length is infinite, the dynamical critical exponent z for CSD is defined 
by the relation 

r oc U . (3) 

The results are shown in figs. 1 and 2. Fig. 1 suggests that at K r eq. (|J) is well 
fulfilled even for small L. The best fit of our data gives 

r(K r ) = 1.20 ± 0.02 . (4) 

At K r /2, fig. 2 suggests that the regime of eq. (|3|) was not yet reached, and that by 
increasing L further the fitted value of r would decrease. We conclude that 

r{K r /2) < 0.79 ±0.09 , (5) 

the number on the r.h.s. being the fitted value for 32 < L < 128. Thus the algorithm 
does indeed perform better at higher temperatures. 

In [II] we analyzed the VMR algorithm for the Discrete Gaussian model. The 
BCSOS algorithm presented here performs similarly to the H-algorithm discussed 
there. Because of the constraints, we did not find a way to include the analogue of 
the I-algorithm that further improved the performance in the case of the Discrete 
Gaussian model. 
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Figure 1: Log-log plot of r vs. L at K = K r . The solid line is the fit with eq. (|3p 
for 8 < L < 128. 



Figure 2: Log-log plot of r vs. L at K = K r /2. The solid line is the fit with eq. @ 
for 32 < L < 128. 



Table 1: The exponent z for fixed L/£ in the smooth phase. 
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Finite size scaling of the autocorrelation time in the smooth 
phase 



In (T2j it is conjectured that the autocorrelation time obeys the usual finite size 
scaling law 

r = L z f(L/0, (6) 

where / is some smooth function. We therefore took data for fixed ratios 
Fortunately, there is an exact formula for the correlation length of the BCSOS 
model (ref. [g], eq. 8.11.24). We made runs for L/£ = 0.25, 0.5, 1, 2, 4 and 8, with 
L chosen to be a power of two, 8 < L < 256. 

For each value of we fitted the values of r, obtained from the autocorrela- 
tions of the one-sublattice energy, with eq. @. In each case the fits were good, but, 
as opposed to eq. @, the value of z was not constant. In table 1 we give the fitted 
values of z as a function of L/£. We see that z slowly increases as L/£ decreases, 
and this increase seems consistent with the value z = 1.20(2) at the transition point 
(L/£ = 0). 

This finding is compatible with a scaling law r oc £ 2 in the limit L — > oo, but 
with an exponent z smaller than 1.20(2). In all cases investigated up to now the 
two exponents were equal, so such a situation would be very interesting. We stress, 
however, that we do not have data close enough to the thermodynamic limit in order 
to claim that the law r oc £ z was well checked. 

In some cases logarithmic corrections to the finite size scaling variable L/C, turn 



out to be important ||13|| . We tried for the same data for r in the smooth phase the 
ansatz 

where we took z = 1.20. We looked for a value of a for which the values of t/L z 
plotted against L/(£(ln£) a ) collapse onto one curve. In fig. 3 we show the result for 
a = —0.7. The range of a for which reasonable collapse occured was —0.9 < a < 
—0.6. Thus, within our precision, this modified finite size scaling is consistent with 
all data and with the value z = 1.20. We caution however that a collapse of the 
data of roughly the same quality can be found for 1.05 < z < 1.25. 

We do not claim that eq. (0) is the "correct" finite size scaling law, but it is 
certainly closer to the truth than eq. (BJ). 
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Figure 3: Finite size scaling plot with z = 1.20 and a = —0.7. Notice the logarithmic 
axes. Different marker symbols are used for different values of L/£. 
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Conclusions and outlook 



We have presented for the first time a successful cluster algorithm for a model with 
constraints. We carefully measured and analyzed the performance of our algorithm. 
We also performed a finite size scaling analysis for the autocorrelation time. 

We would like to point out that a completely different algorithm for vertex 
models was also developed recently |14| . Another research direction we pursue is 
the generalization of the VMR algorithm to more complicated constraints, which 
may be of relevance for comparing to interface experiments. 

While further algorithmic studies are interesting by themselves, we point out 
that the algorithm described in this letter was used in the (to date) most precise 
determination of the roughening (Kosterlitz-Thouless) transition for a variety of 
SOS models (including the dual of the XY model) [p~0f| . 

Finally, we point out that our cluster algorithm can be used, in slightly modified 
form, for accelerating quantum Monte Carlo simulations for XXZ chains [ITS] . 
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